\section{Life insurance with Vasicek interest rate model}
We have a basic two-state life insurance model, where $\kappa$ is receivable upon death before time $T$, and a premium $p$ is paid until the random time $\tau:=$ time of death, or an end time $T$, whichever comes first. Here we take $p=5$, $T=10$ and $\kappa = 1000$.\\

\noindent
Death occurs with intensity

$$
\mu(t) := \lim_{h\rightarrow 0} P\left\{X(t+h)=1|X(t)=0\right\}
$$
where $X(t)$ denotes the state at time $t$. We assume that

$$
\mu(t)=a+bc^{30+t}
$$
where $a=10^{-3}$, $b=10^{-4}$ and $c=1.1$.\\

\noindent
The objective is to evaluate the prospective reserves,

\begin{align}\label{eq:reserve}
V(t):=E\left[R(t)|X(t)=0\right]
\end{align}
where the stochastic process $\{R(t)\}$ is given by

$$
R(t)=\int_t^Te^{-\int_t^sr(u)du}(p1_{\{X(s)=0\}}ds - \kappa dX(s))
$$
Furthermore we assume that the rate of interest is stochastic and governed by the Vasicek model
$$
dr(t)=\alpha(b-r(t))dt+\sigma dW(t)
$$
with $\alpha=0.2$, $\sigma = 0.1$, $b=0.05$ and $r(0)=0.044$.\\

\noindent
Assume that we have simulated the interest rate and want to calculate the reserve at time $t$. Then $r(u)$ can be considered to be deterministic. Recall that $E[1_{\{X(s)=0\}}|X(t)=0] = e^{-\int_t^s\mu(u)du}$ and $E[dX(s)|X(t)=0] = e^{-\int_t^s\mu(u)du}\mu(s)ds$. Then equation~\eqref{eq:reserve} becomes

\begin{align}
V(t)=\int_t^Te^{-\int_t^sr(u)+\mu(u)du}(p - \kappa \mu(s))ds.
\end{align}
which is equivalent to the following differential equation
\begin{align}
dV(t) = ((r(t)+\mu(t))V(t)-p+\mu(t)\kappa)dt,\quad V(T) = 0.
\end{align}

We can now use Euler's (backward) method to evaluate this integral and it results in the following iterative schema

$$
V(t_{i-1})=V(t_i)-((r(t_i)+\mu(t_i))V(t_i)-p+\mu(t_i)\kappa)(t_i-t_{i-1})
$$
and we notice that if we simulate the interest on the same timepoints as we use in the calculation of $V(t)$, then we can easily perform the integration by solving the differential equation.\\

\noindent
Now all we have to do is simulate the interest rate from the Vasicek model. If we do this by Euler's method then it results in the following iterative schema

\begin{align*}
r(t_{i+1})=r(t)+\alpha(b-r(t))(t_{i+1}-t_i)+\sigma \sqrt{t_{i+1}-t_i}Z_{i+1}
\end{align*}
where the $Z_i$'s are N(0,1) random variables that needs to be simulated.\\

\noindent
We have combined this in algorithm \ref{Problem3Algo}, which is implemented in the code in appendix \ref{code3}. We find that the reserve at time zero is approximately 10.0, and to get an understandig of the variability of the estimate we have produced the following histogram, see figure~\vref{fig:reserve}, of $V(0)$'s for 1000 simulated interest rate paths, and one sees that there are some skewness in the $V(0)$'s.

\begin{algorithm}
\caption{Estimation of $V(0)$ with Vasicek interest rate model}
\label{Problem3Algo}
\begin{algorithmic}[1]
  \STATE $V_{sum}\leftarrow 0$
  \FOR {$j \leftarrow 1$ to $N$}
  \STATE $NrSteps \leftarrow 1000$
  \STATE $h\leftarrow\tfrac{10}{NrSteps}$
  \STATE $r_0\leftarrow 0.044$
  \FOR {$i \leftarrow 1$ to $NrSteps$}
  \STATE Generate $Z$ $\sim$ N(0,1)
  \STATE $r_i \leftarrow r_{i-1} + 0.2(0.05-r_{i-1})h+0.1\sqrt{h}Z$
  \ENDFOR
  \STATE $V_{NrSteps}<-0$
  \FOR {$i \leftarrow NrSteps$ to $1$}
  \STATE $V_{i-1}=V_i-((r_i+\mu(ih))V_i-5+\mu(ih)1000)h$
  \ENDFOR
  \STATE $V_{sum}\leftarrow V_{sum}+V_0$
  \ENDFOR
  \RETURN $\tfrac{1}{N}V_{sum}$
\end{algorithmic}
\end{algorithm}
 
\begin{figure}
\centering
\includegraphics[width=10cm]{problem3_1.eps}
\caption{Euler method: Variability in simulated values of the reserve using the Euler scheme for simulating interest rates.}
\label{fig:reserve}
\end{figure}

\noindent
We could also solve this problem in a different way. Instead of moving the expectation inside the integral above we could simulate the time of death as well as the interest rate. Then for each death we can calculate the amount of premium they paid by the insured, and if the death is before $t=T$, then we must subtract the death sum paid to the insured.\\

\noindent
Since we know that the death times have the distribution function $F(x) = 1-e^{-\int_0^x\mu(t)dt}$ we can generate samples from this process by using the inverse transform method. That is we generate $u\sim\textrm{ U(0,1)}$, set $u = 1-e^{-\int_0^x\mu(t)dt}$ and solve for $x$. Since $\tilde{U}=1-U$ where $U\sim\textrm{ U(0,1)}$ is again U(0,1) we can write the equation that we want to solve as $-\log(u) = \int_0^x\mu(t)dt$, which can then be solved numerically for $x$. In our case we have $\int_0^xa+bc^{30+t}dt = ax+c^{30}b\tfrac{c^x-1}{\log(c)}$.\\

\noindent
Given a simulated interest rate we can, once we have simulated a number of death times, calculate the reserve at time 0 for each death time, sum them, and divide the total with the number of simulated death times. This gives an estimate of $V(0)$ for the given simulated interest rate. This can then be performed for a number of simulated interest rates, and the mean of these $V(0)$'s is then an estimate of the true $V(0)$.\\

\noindent
In algorithm~\ref{Problem3bAlgo} pseudocode for performing these operations is given. In this $ValueOfDeathTime$ is the function that calculates the value of a given death time and interest rate. The algorithm as well as the function is implemented in appendix~\ref{code3b}. Upon inspection of the histogram in figure~\vref{fig:reserve_Brute} of the simulated $V(0)$'s, it seems that this method also has some skewness. Furthermore it takes a lot longer to run for fewer simulated interest rate paths since we need to simulate 10000 death times in order to have a resonable estimate of $V(0)$. The reserve estimated in this way is 9.2, and as the histogram shows this methods has considerably larger variance.\\

\begin{algorithm}
\caption{Estimation of $V(0)$ with Vasicek interest rate model}
\label{Problem3bAlgo}
\begin{algorithmic}[1]
  \STATE $V_{sum}\leftarrow 0$
  \FOR {$j \leftarrow 1$ to $N$}
  \STATE $NrSteps \leftarrow 1000$
  \STATE $h\leftarrow\tfrac{10}{NrSteps}$
  \STATE $r_0\leftarrow 0.044$
  \FOR {$i \leftarrow 1$ to $NrSteps$}
  \STATE Generate $Z$ $\sim$ N(0,1)
  \STATE $r_i \leftarrow r_{i-1} + 0.2(0.05-r_{i-1})h+0.1\sqrt{h}Z$
  \ENDFOR
  \STATE $V_{total}<-0$
  \FOR {$i \leftarrow 1$ to $n$}
  \STATE Generate $U$ $\sim$ U(0,1)
  \STATE $deathTime \leftarrow$ find $\left\{x\ge 0|-log(U) = ax+c^{30}b\tfrac{c^x-1}{\log(c)}\right\}$
  \STATE $V_{total} \leftarrow V_{total} + ValueOfDeathTime(deathTime, r)$
  \ENDFOR
  \STATE $V_{sum}\leftarrow V_{sum}+V_{total}/n$
  \ENDFOR
  \RETURN $\tfrac{1}{N}V_{sum}$
\end{algorithmic}
\end{algorithm}

\begin{figure}
\centering
\includegraphics[width=10cm]{problem3_2.eps}
\caption{Histogram of simulated $t=0$ reserves when simlulating both interest rates and time of death.}
\label{fig:reserve_Brute}
\end{figure}

\noindent
We could probably improve upon this method by using that we don't need the detailed information for the insured that survives until time $T$. This means that we could use stratified sampling with two strata; death time larger than $T$ and death time smaller than or equal to $T$. Then we would not need to simulate the death times for the first stratum since the reserve is independent of the actual time of death, and we would only need to estimate the probability of the second stratum. We could then focus on simulating death times smaller than or equal to $T$. This method will however still be slower than the first method, where we used our knowledge of the distribution of the death times to calculate the reserve, especially if we want to estimate the probability of the strata by simulation.

\newpage
